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SUMMARY 


The  application  of  an  integrO'>dif fcrential  approach  in  the  numerical  study  of  unsteady  viscous  flows 
about  airfoils  is  described.  Two  different  procedures  are  presented.  A procedure  based  on  a stream 
function*vorticity  formulation  and  on  a transformation  technique  is  used  in  a study  of  a flow  about  an 
impulsively  started  9X  thick  Joukowski  airfoil  at  an  angle  of  attack  of  15"  and  a Reynolds  number  of  1000. 
Numerical  results  are  presented  and  compared  with  available  finite-difference  results.  A second  procedure 
based  on  a velocity-vorticity  formulation  and  on  a hybrid  finite  difference-finite  element  technique  is  used 
in  a study  of  a flow  about  an  oscillating  12Z  thick  Joukowski  airfoil  at  a Reynolds  number  of  1000.  With 
either  procedure,  the  unique  ability  of  the  integro-dif ferential  approach  to  confine  the  solution  field  to 
the  vortical  region  of  the  flow  is  utilised.  It  is  shown  that  this  ability  offers  great  computational 
advantages . 

1.  INTRODUCTION 

During  the  past  two  decades,  there  has  been  a gradual  but  persistent  increase  in  the  willingness  of 
aerodynamicists  to  accept  the  computer  as  a valuable,  perhaps  indispensable,  tool  in  the  prediction  of  fluid 
flows.  This  willingness  is  a consequence  of  the  impressive  progress  made  in  recent  years  in  the  routine 
numerical  solutions  of  differential  equations  describing  many  types  of  potential  and  boundary  layer  flows. 
For  general  viscous  flows,  i.e.,  flows  past  finite  solid  bodies  and  involving  appreciable  regions  of 
separation,  a capability  for  the  routine  computation  of  the  flows  has  yet  to  be  established,  and  the 
willingness  to  accept  computer  solutions  has  spurred  eatensive  research  efforts  in  this  regard. 

Because  of  its  fundamental  nature  and  practical  importance,  the  general  viscous  flow  problem  has  been 
a focal  point  of  fluid  dynamic  research  for  more  than  a century.  The  mathematical  difficulties  attendant 
to  the  analytical  solution  of  the  equations  describing  the  general  viscous  flow  are  known  to  be  formidable. 
Numerical  methods  offer  at  the  present  the  only  prospect  of  accurste  quantitative  solutions  under  reasonably 
general  circumstances. 

Until  relatively  recently,  computational  aerodynamicists  have  emphasised  the  application  of  the 
f inite-difference  approach  to  the  general  viscous  flow  problem.  This  emphasis  is  an  outgrowth  of  the  success 
enjoyed  by  the  finite-difference  tnethods  in  the  solution  of  various  boundary  layer  problesu.  Although  the 
boundary  layer  equations  are  simpler  than  the  Navier-Stokes  equations  describing  the  general  viscous  flow, 
the  two  types  of  equations  are  both  non-linear  partial  difference  equations.  Furthermore,  the  asMaentum 
equation  for  the  time-dependent  general  viscous  flow  is  parabolic,  as  is  the  boundary  layer  momentum 
equation.  It  had  been  hoped,  therefore,  that  the  finite-difference  approach  would  lead  to  a routine 
computational  capability  for  the  time-dependent  general  viscous  flow,  with  steady  flow  solutions,  when  they 
eaist,  considered  as  asymptotic  limiting  solutions  at  large  time  levels.  In  reality,  however,  the 
application  of  the  finite-difference  approach  to  the  genera]  viscous  flow  problem  has  fallen  far  short  of 
earlier  expectations.  Host  of  the  general  viscous  flows  of  practical  importance  in  aerodynamics  today  remain 
beyond  the  scope  of  the  finite-difference  approach. 

Several  years  ago,  the  first  author  of  this  article  identified  several  major  obstacles  that  caused  the 
disparity  of  successes  oi  the  finite-difference  approach  as  applied  to  the  boundary  layer  and  to  the  general 
incosN)ressible  viscous  flows.  Subsequently,  a program  of  research  was  initiated  with  the  goal  of  developing 
new  approaches  transcending  the  limitations  of  the  finite-difference  approach.  This  research  program  has 
noM  reached  a reasonable  stage  of  completeness.  Progress  made  to  date  in  applying  a new  approach  to  the 
stody  of  unsteady  incompressible  laminar  flows  is  sumarised  in  this  paper.  Selected  numerical  results  are 
presented  for  flows  past  airfoils  to  illustrate  the  application  of  this  new  approach.  Additional  numerical 
results  as  well  as  detailed  descriptions  of  the  specific  nusMrical  procedures  are  presented  in  doctoral 
dissertations  of  the  second  author  (Ref.  1)  and  of  the  third  author,  the  latter  still  being  prepared.  A 
major  emphasis  of  the  present  article  is  to  bring  into  focus  the  authors'  understanding,  acquired  during  the 
past  few  years  of  research,  of  the  interplay  between  the  numerical  and  physical  aspects  of  the  general 
viscous  flow  problem. 

The  distinguishing  feature  of  the  present  approach  is  the  formulation  of  the  general  viscous  flow 
problem  as  a set  of  integro-differential  equations.  This  feature  represents  a major  departure  from  finite- 
difference  methods  based  on  the  formulation  of  the  problem  as  differential  aquations.  Various  attributes  of 
the  integro-differential  approach  have  b«Mn  studied  and  described  in  earlier  articles  (Refs.  1 to  11).  In 
particular,  it  has  been  conclusively  demonstrated  by  analyses  and  numerical  illustrations  that  the  approach 
possesses  a unique  ability  to  confine  the  solution  field  to  the  vortical  region  of  the  flew.  For  moat 
problems  of  practical  interest,  the  flow  Reynolds  number  is  high  and  the  vortical  region  represents  only  a 
ssMll  fraction  of  the  total  ftowfiald.  The  ability  to  confine  the  solution  field  to  the  vortical  region 
therefore  offers  the  possibility  of  a drastically  i^roved  computational  efficiency. 

Although  the  present  paper  is  concerned  only  with  unsteady  incompressible  laminar  flows,  the  utility 
of  the  integro-differential  approach  is  not  restricted  to  such  flows.  In  a recent  article  (Ref.  3),  the  use 
of  the  integro-differential  approach  in  conjunction  with  a two-equation  turbulence  modal  was  demonstrated. 
Turbulent  flow  results  obtained  were  found  to  agree  well  with  experimental  data.  An  integro-differential 
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formulation  was  presented  for  compressible  flows  three  years  ago  (let.  4).  Because  of  the  limited  resources 
availablei  however,  efforts  of  the  present  authors  in  implementing  the  integro-dif ferential  approach  have 
been  limited  to  incompressible  flow  problems.  It  has  also  been  found  that  the  general  viscous  flow  problem 
may  be  formulated  as  a set  of  integral  representations  (Ref.  5)  of  field  variables.  The  solution  procedures 
developed  for  the  integro-differential  approach  are  directly  applicable  to  methods  based  on  the  integral* 
representation  approach.  In  particular,  for  steady  flows,  a method  based  on  the  integral  repraaentetion 
approach  has  been  established  (Ref.  6)  and  a general  computer  code  for  internal  flow  problems  is  being 
prepared. 

Recent  literature  contains  several  articles  by  other  researchers  describing  results  obtained  using 
the  integro-differential  approach  (Refs.  12  to  15).  The  relatively  large  amounts  of  computing  needs 
associated  with  these  efforts  as  well  as  with  our  earlier  work  (Ref.  7)  led  to  the  opinion  that,  even  though 
the  integro-differential  approach  permits  the  solution  field  to  be  confined  to  the  vortical  ragiem  of  Che 
flow,  the  approach  is  computationally  inefficient  (Ref.  16).  This  impression  is  misleading.  The  present 
authors  have,  in  the  past  few  years,  progressed  through  several  stages  of  development  of  the  iekegro- 
differential  approach;  each  stage  has  led  to  a substantial  reduction  in  computing  need.  The  present 
computing  need  is  about  1/20  of  that  reported  in  Raf.  7 and  for  two-dimensional  problems  is  about  1/5  of  that 
required  by  Che  more  efficient  finite-difference  methods  available  today.  The  reduction  in  needed  computer 
time  is  expected  to  be  much  greater  for  three-dimensional  problems. 

Two  different  numerical  procedures  are  outlined  in  this  paper.  A procedure  uciliting  a transforma- 
tion technique  is  used  in  a study  of  a flow  about  an  impulsively  started  9X  thick  airfoil.  A second  procedure 
utilizing  a hybrid  finite  difference-finite  element  technique  is  used  in  s study  of  flows  about  an 
oscillating  12Z  chick  airfoil. 

2.  KINETICS  AND  KINEMATICS  OF  THE  PLOW 

The  time-dependent  Navier-Stokes  and  continuity  equations  for  a fluid  with  constant  density  p , 
constant  kinematic  viscosity  v , and  subject  to  negligible  body  forces  are  well-known  and  are  expressible  in 
terms  of  the  velocity  v and  pressure  p as: 

♦ (V  .7)  ; - 


\>V*v 


(1) 

(2) 


where  t is  the  time. 

Equations  (1)  and  (2)  are  in  principle  sufficient  for  Che  determination  of  v and  p,  provided  that 
adequate  initial  and  boundary  condicions^are  known.  The  boundary  conditions  occuring  most  frequently  are  the 
"no-slip"  condition  at  t^e  surface  of  solid  bodies  over  which  the  fluid  passes.  For  problems  involving  an 
infinite  fluid  domain,  v and  p at  infinity  must  also  be  specified. 

It  is  convenient  to  introduce  the  vorticity  vector  u)  defined  by 

5 » V - I (3) 

and  to  consider  the  "vorticity  transport"  equation 

|S  - ^*(»  , 2)  ♦ (4) 

obtained  by  taking  the  curl  of  both  sides  of  Eq.  (1)  and  using  Bqa.  (2)  and  (3). 

The  sec  of  equations  (2)^to  (4),  with  v and  (i>  as  dependent  variables,  replaces  the  eq^t  of 
equations  (1)  and  (2)  in  idtich  v and  p arc  dependent  variables.  There  are  several  reas^s  for  using  u in 
the  formulation  of  the  problem.  A principal  reason  is  that  the  formulation  in  terms  of  (s  allows  clear  and 
separate  delineation  of  the  kinetic  and  kinematic  aspects  of  the  problem.  As  a consequence  of  this  clear 
delineation  several  important  characteristics  of  the  general  viscous  flow  problem,  not  obvious  in  the  p and 
v formulation,  become  menifest.  Major  obstacles  to  the  finite-difference  solution  of  the  general  viscous 
flow  problem  are  then  traceable  to  physical  processes  involved  in  the  development  of  the  flow  patterns.  The 
identification  of  the  physical  origin  of  these  obstacles  is  of  course  a prerequisite  for  the  establishment  of 
solution  methods  that  overcome  these  obstacles. 

♦ 

The  kinetic  e.pect  of  the  problea  deele  with  the  chenge  of  the  vorticity  field  u with  tiae.  Thie 
especC  ie  deecribed  by  Eg.  (4).  For  an  inviecid  fluid,  the  leet ^ora  in  Iq.  (4)  veniihee  nd  the  vorticity 
ie  convected  with  the  fluid  in  the  eenee  that  the  vorticity  flux  in  .nde,  tdiere  n ie  a unit  noraal  vector  of 
the  aurface  eleaent  da,  aaaociated  with  each  aatarial  aleaant  da  aoving  with  the  fluid  reaaina  a cooatant 
for  all  tiaea.  Thia  well-known  thcorea  of  Halaholca  ia  aodified  in  tha  caaa  of  a real  fluid  by  the  proceta  of 
vorticity  diffuaion.  According  to  Eq.  (4),  changea  in  tha  flux  of  vorticity  acroaa  a aurface  alaaant  that  ia 
aoving  with  the  fluid  and  ia  in  the  interior  of  the  fluid  doaain  takea  place  only  by  diffuaion.  Vorticity 
flux  cannot  be  created  in  the  interior  of  a fluid. 

If  the  fluid  ia  in  contact  with  aolid  bodiaa  in  notion,  the  no-alip  condition  providea  a aaehaniaa  for 
tha  generation  of  vorticity  at  tha  aolid  aurfacea.  In  the  caaa  where  tha  fluid  it  inUially  at  reat  and 
occupiea  an  infinite  dnuin  bounded  internally  by  aurfacea  of  aolid  bodiaa  alao  initially  at  reat,  the 
vorticity  ia  obvioualy  aero  avtryidiere  initially,  lubaequant  aotion  of  the  aolid  bodiaa  aeta  up  corratpond- 
iog  notion  of  the  fluid.  laadiately  after  tha  oaaat  of  tha  notioa,  ainca  tha  vorticity  flux  in  the  interior 
of  the  fluid  doaain  can  change  only  aa  a contaquaace  of  diffuaion,  the  vorticity  ia  everywhere  aero  except  at 
the  aolid  aurfacea.  That  ia,  tha  flew  iaadiataly  after  tha  oaaat  of  the  ■etiea  haa  a non-aaro  tangential 
velocity  relative  to  the  aolid  bodiaa  at  the  aolid  hoaadariea . Tha  diacontiaaity  in  tangential  velocity 
rapraaenta  a aheet  of  concentrated  vorticity  at  the  aolid  houndariaa.  At  anbaaquant  tiaM  lavnla,  thia 
concentrated  vorticity  at  tha  aolid  beundariaa  apraada  into  tha  interior  of  tha  fluid  doaain  by  diffuaion 
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•nd,  once  there,  !•  transported  avay  froa  the  boundaries  by  both  convection  and  diffusion.  Ac  the  eaae  tiae, 
vorticiCy  is  continually  generated  on  the  solid  boundaries.  The  general  flow  pattern  therefore  includes 
vortical  regions  surrounding  the  solid  bodies  and  vortical  wake  regions  trailing  the  solid  bodies.  Outside 
of  these  vortical  regions,  the  flow  is  essentially  free  of  vorticiCy  and  is  therefore  potential.  If  Che 
Reynolds  nuaber  of  Che  flow  is  not  small,  then  the  vorticiCy  spreads  by  diffusion  only  a short  distance  froa 
the  solid  surface  before  being  carried  away  with  the  fluid.  Thus  there  exists  a large  region  of  the  fluid, 
ahead  and  to  the  side  of  the  solid  bodies,  which  is  free  of  vorticiCy.  In  fact,  at  any  finite  tiae  level 
after  the  onset  of  the  motion,  the  vortical  regions  are  finite  in  extent  while  the  potential  flow  region  is 
infinite  in  extent. 

The  kinematic  aspect  of  the  problem  relates  the  velocity  field  at  any  given  instant  of  tiae  Co  the 
vorticity  field  at  chat  instant.  This  aspect  is  described  by  Eqs.  (2)  and  (3)  together  with  appropriate 
velocity  boundary  conditions  on  the  solid  surfaces  and  at  infinity.  If  the  velocity  field  is  known  at  any 
given  inatant  of  time,  Che  corresponding  vorticity  field  is  iasediately  obtainable  by  differention,  using 
Eq.  (3).  If  the  vorticity  field  is  known,  Che  velocity  field  is  obtainable  by  integrating  Bq.  (3),  subject 
Co  the  solenoidal  condition,  Eq.  (2).  One  method  of  evaluating  % is  to  Cake  Che  curl  of  Eq.  (3)  and,  by 
using  Eq.  (2),  obtain  a vector  Poisson's  equation  for  v in  Che  form 

V*v  “ X I (5) 

Alternatively,  Eq.  (2)  implies  the  existence  of  a vector  potential  $ which  is  related  to  v by 

^ ^ X ? (6) 

Equation  (3)  and  (6)  gives  s vector  Poisson's  equation: 

V X V X ? - u (7) 

A general  solution  procedure  based  on  the  kinetic**kineMtic  formulation  which  enables  the  pattern 
flow  development  to  be  followed  computationally  is  as  follows.  With  known  spatial  distributions  of  v and^U) 
at  a given  time  level,  the  kinetic  aspect,  i.e.,  Eq.  (4)  together  with  appropriate  boundary  conditions  on  u) , 
is  Created  numerically  and  a new  distribution  of  vorticity  at  a subsequent  time  level  is  established.  The 
kineMtic  aspect,  i.e.,  Eqs.  (2)  and  (3),  or  Eq.  (7),  together  with  appropriate  boundary  conditions  for  v, 
or  T , is  Chen  Created  and  a new  velocity  distribution  corresponding  to  the  new  vorticity  distribution  at  the 
new  time  level  is  obtained.  This  completes  a computational  loop  to  advance  the  solution  from  an  old  time 
level  Co  a new  one.  Repeated  applications  of  this  computational  loop  yields  a transient  solution.  In  the 
limit  of  large  time,  if  the  velocity  boundary  conditions  are  independent  of  tisM,  Chen  either  a steady  state 
solution  or  a periodic  vortex  shedding  solution  is  obtainable. 

The  implementation  of  Che  above  outlined  solution  procedure  consists  of  a number  of  component 
problems  as  shown  in  Figure  1.  The  literature  contains  a variety  of  numerical  methods  Chat  differ  from  one 
another  in  the  specific  techniques  for  treating  the  component  problems.  The  general  framework  for  these 
numerical  jMthods,  however,  are  identical.  ^If  the  "primitive”  variable  p is  used  in  place  of  the  "derived" 
variable  u>»  or  if  y is  used  in  place  of  v,  the  component  problems  are  then  concerned  with  the  computation 
of  these  different  variables.  The  spirit  of  the  general  soluCion-procedure  framework,  however,  remains  Che 
same.  The  major  difficulties  encountered  by  researchers  can  be  conveniently  reviewed  by  a critical 
examination  of  the  component  problems. 

3.  MAJOR  OBSTACLES  AND  ALTERNATIVES 

In  this  section,  several  major  obstacles  experienced  in  previous  finite*dif ference  solutions  of  the 
general  viscous  flow  problem  are  reviewed.  As  stated  earlier,  the  finiCe-dif fcrence  approach  enjoys 
remarkable  success  in  the  prediction  of  boundary  layer  flows.  The  present  obstacles  are  traceable  Co  the 
important  differences  between  the  boundary  layer  flow  and  the  general  viscous  flow  in  flow  patterns  and  in 
Che  physical  processes  of  flow  development.  It  is  worthy  of  note  chat  except  the  establishment  of  the 
Initial  solution,  each  of  the  component  probleim  for  the  general  viscous  flow,  as  shown  in  Fig.  I,  present 
difficulties  that  are  absent  in  the  boundary  layer  flow.  These  difficulties  dictate  the  development  of 
innovative  approaches  for  the  general  viscous  flow  problem.  Several  such  approaches  are  described  here. 

3.1.  Grid  System-'-Coordinate-Transformation  and  Pinite^Element  Methods 

As  discussed  earlier,  the  flowficld  in  general  consists  of  vortical  regions  surrounded  by  a potential 
flow  region.  For  high  Reynolds  number  flows,  if  no  appreciable  flow  reversal  occurs  near  Che  solid  bodies, 
the  vorticity  spreads  to  Che  sides  of  the  solid  bodies  mainly  by  diffusion.  Since  the  diffusive  process  is 
relatively  slow  at  high  Reynolds  numbers,  Che  vortical  regions  near  the  solid  bodies  are  chin  and  Che 
boundary  layer  simplifications  are  justifiable  there.  The  vortical  wakes  derive  their  vorticity  from  the 
Chin  boundary  layers  and  are  therefore  also  chin.  The  success  of  the  boundary  layer  theory  is  attributable 
to  Che  fact  chat  the  effects  of  these  vortical  layers  on  the  potential  flow  may  be  treated  as  a relatively 
small  perturbation  to  the  potential  flow  past  the  solid  bodies  in  the  absence  of  the  vortical  layers.  This 
fact  permits  the  potential  region  Co  be  studied  separately  from  the  vortical  regions.  In  particular,  it 
becomes  possible  to  ascabliah  the  flow  velocity  immediately  adjacent  to  the  boundary  layer  through  potential 
flow  analyses,  employing  iterative  or  smtching  procedures  if  needed  to  properly  account  for  the  effects  of 
Che  vortical  layers.  Consequently,  it  is  possible  Co  confine  the  solution  field  of  the  boundary  layer 
equation  Co  the  thin  vortical  layers  near  the  solid  bodies. 

In  general  viscous  flows,  appreciable  regions  of  flow  reversal  occur  near  the  solid  bodies.  The 
vorticiCy  now  spreads  to  the  sides  of  the  solid  bodies  by  both  diffusion  and  convection.  The  vortical 
regions  are  not  thin  and  their  presence  "strongly"  influences  the  potential  flow.  With  the  finite-dlfference 
approach,  the  solution  field  must  include  the  potential  region  as  well  as  the  vortical  regions.  Matching 
procedures  are  not  eaaily  applicable. 

For  attached  boundary  layers,  it  is  permissible  to  employ  an  orthogonal  coordinate  ayatem  with  one  of 
Che  coordinate  sms  coinciding  with  the  solid  surface,  fuch  a coordinate  system  is  curvilinear.  However, 
since  the  solution  field  is  restricted  to  the  thin  region  of  the  boundary  layer,  the  scale  factor  of  the 
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curviline«r  tytccm  ch«ng««  n€gligibly  scrotf  Che  eoluCion  field.  The  differeoCiel  equeciooe  deecribing  Che 
flow,  ^en  expreeted  in  Ceras  of  Che  curvilinear  coordineCee,  ere  foraelly  identicel  with  Che  equecione 
expreaeed  in  Che  Certeeien  coordineCee,  che  curveCure  cenw  being  negligibly  eaell  within  Che  contest  of 
boundery-leyer  epproxiaecioae . Ae  e coneequence,  it  ie  eCreighC forward  to  use  e finite*dif ference  grid 
■yeten  with  the  airfoil  eurfece  represented  by  e grid  line,  luch  e grid  eystea  ie  "boundery-fiCCed"  in  the 
sense  Chet  date  points  on  the  solid  boundaries  ere  pervitted  to  coincide  iHch  grid  points. 

For  the  general  viscous  flow  problea,  che  solution  field  is  not  thin  end  the  scale  factor  in  a 
curvilinear  coordinate  systea  fitted  to  che  solid  boundaries  can  change  substantially  across  the  solution 
field.  The  differential  equations  describing  the  flow,  when  aspressed  in  the  curvilinear  coordinates, 
contain  curvature  teras  that  are  not  netUgible.  The  use  of  a boundary •fitted  finite^dif farance  grid  is 
considerably  acre  coaplas  than  it  is  in  che  case  of  an  attached  flew. 

Much  of  che  earlier  works  using  the  finite-difference  approach  to  treat  the  general  viscous  flow 
problem  atceapCed  to  avoid  the  coaplexity  introduced  by  Che  use  of  curvilinear  coordinata  systea.  In  these 
works,  the  differential  equations  are  expressed  in  Cartesian  Coordinatas  and  a rectangular  grid  systea  in  Che 
physical  plana  is  used.  Such  a systea  is  not  boundary*fitCed.  That  is,  grid  points  representing  the 
boundary  in  general  do  not  coincide  with  che  solid  boundaries,  and  intarpolation  procedures  are  needed  to 
accoaaodate  Che  boundary  geoaetry. 

Various  authors  treating  different  types  of  flow  problems  (tefs.  17  to  21)  ere  in  agreement  that  the 
proper  handling  of  the  boundary  conditions  are  of  dominant  laportance  in  Cha  numerical  solution  of  flow 
probleas.  existing  evidence  points  to  the  use  of  grid  systems  that  are  not  bodyfitted  as  one  of  cha 
culprits  for  the  very  liaiCed  capability  of  earlier  aethods  in  treating  the  general  viscous  flow  problea.  As 
is  wall-known,  the  Naviar-Stokes  equations  describing  Che  general  viscous  flow  are  coMon  to  an  astonish- 
ingly rich  variety  of  flow  phooomono  drosticoLlydifforin,  frow  ono  onochornot  only  qunnticotiwoly  but  also  in 
charactor.  Tha  draatic  diffarancaa  ara  raaulta  of  diffaraneaa  In  tka  (aoaacriai  and  tha  iayoaad  condltlona 
at  tha  flow  boundariaa.  Claarly,  aecurata  rapraaantatlon  of  tha  flow  boundarlaa  and  of  tha  boundary 
condition!  ara  of  paraaount  iaportanca  in  tha  pradictlon  of  ,anaral  viacoua  flows. 

Boundary- fittad  grid  aystaaw  ara  obtainabla  by  tha  uaa  of  coordinate  transfoTaation  aathoda. 
hltarnativaly,  tha  finita-alaaant  aathoda  (Baf.  IB)  accoapliah  tha  aaaa  purpoaa.  With  aithar  approach,  tha 
ralativa  aiaplicity  of  tha  Cartaaian  finita-dif faranca  raprasantation  ara  bartarad  for  a battar  nuaarical 
raprasantation  of  tha  boundary  gaoaatry. 

Until  racantly,  auccaaaful  nuawrical  aolutiona  of  tha  ganaral  viacoua  flow  probloa  hava  baan  obtained 
only  for  those  boundary  gaoaatrias  for  which  convenient  analytical  tranaforaatlon  foraulaa  ara  available. 
In  perticular,  conformal  trana format  ions  which  permit  straight  forward  calculation  of  tha  scale  factor  for 
tha  curvilinear  coordinata  syatam  and  of  the  curvature  terms  in  the  differential  equationa  hava  bean  utilised 
by  aavaral  rasaarchara.  In  recant  years , numerical  mathods  for  tha  astabliahmant  of  boundary-fitted 
coordinate  syatema  have  been  developed  (Raf.  17).  The  numarleal  transformations  need  not  ba  conformal  and  is 
not  rastrictiva  regarding  tha  boundary  gaomatry. 

Tha  baa.c  finita-alaaant  aathoda  involves  tha  following  coaponants:  (a)  tha  raforaulatian  of  the 
diffareacial  aquations  as  integral  ralatloaa,  (b)  Cha  division  of  cha  solution  fiald,  l.a.,  cha  ragion  of 
intagratism,  imco  small  shbragions,  called  alamencs,  af  arbitrary  sisaa  and  shapes,  (e)  Cha  usa  of 
iatarpolaclon  functions,  usually  polyncaials,  to  anprass  cha  dapamdamC  variablas  in  each  alamaac  la  Caras  of 
chair  values  at  sslactsd  data  points,  callad  nodas,  aasociatad  with  tha  alaaant,  and  (d)  Cha  evaluation  of 
cha  intagrals  over  aacb  slasHnt,  yialdiag  a aat  of  algebraic  aquations  comtaiaiag  tha  nodal  values  of  tha 
dapandant  variablas  as  umkaouns.  tinea  tha  alamaats  ara  af  arbitrary  sisas  and  shapas,  it  is  straightforward 
to  davisa  a boundary-fictad  syatam  of  elamants.  Tha  raaultiag  sat  of  algebraic  aquations,  howavar,  poasassaa 
coefficient  macrieas  more  ccmplax  than  choss  of  finita-dif fsrsaca  aquations. 

In  tha  present  work,  both  cha  transformation  method  and  cha  finita-alananc  mathod  ara  utlliasd  in 
conjunction  with  Che  Incagro-dlffarantial  approach  in  studies  of  viscous  flows  about  airfoils.  Tha  atudiaa 
danonscratsd  Chat  both  aathoda  ara  sffacclva  in  raaoviag  tha  difficulty  of  accurate  nuaarical  raprasantation 
of  cha  solid  boundaries.  Tha  added  algobralc  conplsaity  of  aithar  method  is  not  aacasaiva  for  two- 
diaanaioaal  flows.  Thera  asist,  howavar,  aavaral  other  dlfficultias  that  ara  not  ramovad  by  aithar  method. 
Thasa  dlfficultias  ara  daacribed  below. 

3.2.  Torticit2_8alution  - Klnatics 

a major  sourca  of  difficulty  in  cha  kinacie  part  of  ganaral  viscous  flow  computation  is  tha 
simulcanaoua  occurrence  of  tha  diffusiva  and  convsetiva  procassas  in  tha  flow.  Those  two  procaasaa  proceed  at 
drastically  different  rates  at  hl|^  Reynolds  aumbars.  Tha  laagth  and  tlma  scales  of  thasa  procaasaa  ara 
vastly  different.  This  fact  loads  to  vastly  diffsrant,  ofcsa  overly  stringent,  raquiramanta  on  tha  tins 
Incarval  and  grid  spacing  bacausa  of  accuracy  and  stability  considaracions.  In  comparison  to  problams  in 
idiich  aithar  cha  diffusiva  procass  or  tha  convactivn  process  is  alone  important,  Cha  kinatic  aspact  of  cha 
ganaral  viscous  flow  is  Immsnsaly  mors  complaa.  far  boundary  layer  flows  in  cwo-dimnaioaa,  cha  diffusiva 
procass  is  of  nagligibla  isvortanca  in  comparison  to  the  coavactiva  procass  along  tha  diraetion  parallel  Co 
cha  body  surface.  In  tha  direction  parpandieular  to  tha  body  surface,  cha  convactlvs  and  diffuaivn  procassas 
ara  of  eoaparsbla  importanca  and  tha  length  ssalas  af  tha  two  procassas  ara  not  vastly  diffaranc. 
Consaquancly,  tha  boundary  layer  aomantun  aquatiams  do  not  soataln  che  above  dascrlbod  diffisulty  asaociatad 
with  cha  ganaral  viscous  flow. 

In  Raf.  3,  an  integral  rapresamtatiao  for  tha  vartiaity  vactar  is  presontod.  This  intagral 
raprasaatacioa  is  eomplataly  aquivalanc  to  tha  dlffarontial vmrtlcity  transport  aquation  (4).  A aolutiam 
pracadura  based  on  this  mprasentatloa  la  suggaatad  U Baf.  S.  This  proaadura  has  baan  callbratad  racantly 
using  simpla  problasa  with  km sum  amalytiaal  nalntlaaa.  A natawarchy  fnatura  af  this  praeadmro  la  that  tha 
kiaacia  procassas  af  diffaslan  and  aonvoctian  ara  raprasamtad  by  aaparata  Uugrals  and  tha  aaatrlbntiama  af 
thasa  ewa  pracassas  to  tha  tlma  varlatian  af  vartiaity  naa  ba  avalaatad  aaparatnly  ia  a caavaaiaat  maaaar. 
This  aaw  praaadurs  is  baiag  davalapad  far  aamplaa  tima  dapandant  problasm.  Tha  work  rapartad  ia  tha  prnsant 
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paper,  however,  are  baaed  on  extenaiona  of  exiating  finite-difference  and  finite-element  methoda.  For  the 
impulaively  started  airfoil,  a method  of  flowfield  segmentation  is  developed  and  used  in  conjunction  with  a 
modified  strongly  implicit  procedure.  For  the  oscillating  airfoil,  a hybrid  finite  difference-finite 
element  method  ia  eatabliahed.  These  mechods  will  be  outlined  in  later  sections  of  this  paper. 

3.3.  Kinematics  - Integral  Representation  and  Flowfield  Segmentation 

It  is  well  known  that  the  computer  time  requirement  for  the  finite-difference  solution  of  a general 
viscous  flow  problem  in  general  increases  rapidly  with  the  flow  Reynolds  number.  The  basic  cause  of  this 
difficulty  is  that  the  kinematic  aspect  of  the  general  viscous  flow  problem  is  described  by  elliptic 
differential  equations,  e.g.  Eq.  (5)  oj  Eq.  (7).  The  finite-difference  solution  of  v at  any  given  point  in 
the  flowfield  ^ependa  on  the  value  of  v at  the  neighboring  points.  It  is  therefore  not  possible  to  evaluate 
the  value  of  v explicitly,  point  by  point,  using  a finite-difference  Mthod.  The  solution  field,  in  fact, 
should  be  infinite  so  as  to  satisfy  the  freestream  condition  at  infinity.  For  numerical  solution,  the 
infinite  solution  field  needs  to  be  represented  by  a finite  number  of  grid  points.  Substantial  inaccuracy 
can  result  if  one  truncates  the  solution  field  and  specifies  freestream  velocity  condition,  or  some  other 
condition,  on  some  outer  boundaries  at  finite  distances  from  the  solid  boundary  (Reference  8).  In  many 
previous  works,  this  outer  boundary  is  placed  at  a distance  of  ten  chord  lengths  or  more  from  the  solid 
boundary  in  hope  that  the  inaccuracy  resulting  from  the  replacement  of  an  infinite  region  by  a finite  one  is 
then  unimportant.  Such  practices  may  lead  to  fundamental  difficulties  and  sxist  be  handled  with  great 
caution.  For  example,  it  can  be  shown  that  if  the  steady  state  equations  are  to  be  treated  and  the  freestream 
velocity  boundary  condition  is  applied  at  a closed  outer  boundary  at  finite  distances  from  the  airfoil,  then 
no  solution  for  an  airfoil  with  non-zero  lift  force  is  possible. 


As  discussed  earlier,  for  high  Reynolds  number  flows,  the  vortical  regions  surroundina  and  trailing 
the  solids  are  relatively  small.  The  effecta  of  viscosity  are  important  only  in  this  region.  The  gradients 
of  field  variables  are  large  and  the  corresponding  length  scale  for  this  region  is  siaall,  particularly  in  the 
boundary  layer.  Consequently,  in  order  to  have  sufficient  solution  resolution  in  this  region,  the  grid  lines 
must  be  very  finely  spaced.  If  this  fine  spacing  is  continued  into  the  potential  region  of  the  flow,  which 
has  been  made  finite  but  is  still  large,  a gigantic  number  of  grid  points  invariably  enter  the  computation 
procedure.  This  implies  that  an  excessive  amount  ofcomputer  time  is  needed  to  solve  the  equations.  With 
increasing  Reynolds  number,  the  length  scale  in  the  vortical  region  decreases.  The  number  of  grid  points, 
and  hence  also  the  needed  computer  tisie,  therefore  increases  with  increasing  Reynolds  number. 

The  coordinate  transformation  and  finite-element  methods  provide  a degree  of  relief  from  the 
excessive  computing  needs  at  high  Reynolds  numbers.  For  example,  Ref.  17  describes  a method  which  maps  the 
fluid  region  into  the  interior  of  a rectangle.  A uniformly  spaced  grid  system  is  used  in  the  transformed 
plane  to  numerically  solve  the  transformed  differential  equations.  In  the  physical  plane,  the  corresponding 
spacing  between  grid  lines  increases  as  the  distance  from  the  solid  boundaries  increases.  Vith  the 
"expanding'*  grid  system  in  the  physical  plane,  it  is  possible  to  substantially  reduce  the  number  of  grid 
points  entering  the  solution  procedure  and  to  make  this  number  less  sensitive  to  the  Reynolds  nuc^er. 

With  a f inite-eleomnt  method,  since  the  eleownts  can  be  of  arbitrary  sizes  and  shapes,  an  expanding 
grid  system  is  not  difficult  to  devise.  To  researchers  accustomed  to  the  finite-difference  approach,  the 
finite-element  procedure  for  setting  up  a set  of  algebraic  equations  may  appear  to  be  circuitous.  It  should 
be  noted,  however,  that  the  coordinate  transfomration  approach,  usually  considered  an  improved  finite- 
difference  procedure,  also  would  appear  to  be  circuitous  to  those  accustomed  to  the  use  of  uniformly  spaced 
finite-difference  grid  in  the  physical  plane.  The  added  complications  of  coordinate  transformation,  or  of 
the  finite-element  approach,  hoivever,  are  amply  justified  by  their  ability  to  accurately  represent  the 
boundary  geometry,  as  discussed  in  Section  3.1.  The  ability  to  have  an  expanding  grid  system  is  an 
additional  important  property  of  these  approaches.  However,  these  approaches  by  themselves  do  not  permit  an 
explicit,  point  by  point,  evaluation  of  the  velocity  values,  and  therefore  the  solution  field  needs  to 
include  the  potential  region  as  well  as  the  vortical  regions  of  the  flow. 

As  discussed  in  Section  3.1,  in  treating  the  attached  flow  problem,  the  solution  field  for  the 
boundary  layer  equations  can  be  confined  to  the  thin  vortical  layers  near  the  solid  bodies.  There  is  no  need 
to  concurrently  solve  the  potential  and  the  vortical  flow  problems.  Rather,  the  t%ro  problems  are  solved 
individually  and  matched.  With  attached  flows,  therefore,  the  difficulty  of  simultaneously  accomendeting 
the  potential  and  vortical  regions  ia  absent.  ThiB  fact  ia  a major  reason  contributing  to  the  success  of  the 
finite-difference  approach  in  solving  boundary  layer  problems. 

For  the  ^neral  viscous  flow  problem,  the  kinetic  aspect  of  the  computation  can  be  confined  to  the 
vortical  regions  of  the  flow.  From  Eq.  (4),  it  is  obvious  that  if  the  vorticity  and  velocity  fields  in  the 
vortical  regions  are  known  at  any  given  time  level,  then  the  time  rate  of  change  of  vorticity  at  that  time 
level  is  determinate.  There  is  no  need  to  go  beyond  the  vortical  regions  and  their  lasediate  neighborhood  to 
establish  a new  vorticity  distribution  at  tha  subsaquant  time  level.  Consequently,  if  the  kinematic  aspect 
of  the  computation  can  also  be  confined  to  the  vortical  region,  then  the  solution  field  need  not  extend  into 
the  potential  flow  region  and  tha  difficulty  of  simultaneously  accoModating  Che  potential  and  the  vortical 
regions  is  removed.  In  Ref.  2,  it  is  shown  that  tha  confinement  of  the  kinematic  aspect  of  the  computetion 
to  the  vortical  region  is  possible  if  the  differential  equations  describing  Che  kinematics  of  the  general 
viscous  flow  problem  is  recast  into  an  integral  representation  of  the  velocity  vector  (or  of  the  vector 
potential),  la  thia  research,  the  integral  repreaentation  for  the  kinematics  of  the  problem  is  used  together 
with  the  differential  vorticity  transport  equation  (4).  The  formulation  of  the  overall  problem  is  named  the 
integro-differeotial  formulecion. 


In  Ref.  2,  it  is  shown  that  Rqs.  (2)  and  (3),  together  with  velocity  boundary  conditions,  can  be 
recast  into  the  following  integral  representation  for  the  velocity  vector: 
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'**'•*’•  r i*  th«  poaition  vector,  u ^ia  th«  vorticlty  vactor,  R ia  any  region  occupied  by  the  fluid  or  the 
aolid,  1 ia  the  boundary  of  I,  and  n ia  the  unit  normal  vector  directed  outward  from  S,  A ■ 47f  and  d • 3 
in  Chree~*dimanaionaI  problema  and  A ■ 2w  and  d • 2 ^in  cvo-dimenaional  problema.  The  aubacripta  ”o"  in 
^K^and  dB^  indicate  that  the  integrationa  are  in  the  r^  apace. 

Bduation  (d)  ia  valid  for  all  poaition  vectora  ia  R and  on  b and  ia  completely  equivalent  to 
Rqa.  (2)  and  (3)  together  with  the  velocity  condition  on  the  boundary  B.  The  uae  of  Rq.  (8)|  in  place  of 
Sqa.  (2)  and  (3),  parmita  the  expliciti  point  by  point,  computation  of  the  velocity  diatribution  in  R, 
provided  that  the  vertieity  diatribution  in  R and  the  velocity  condition  on  B are  known. 


The  aacond  integral  in  Bq.  (8)  repreaenta  the  contribution  of  the  velocity  boundary  condition  to  the 
velocity  field  within  the  boundary.  For  the  problem  of  an  external  flow  peat  a aolid  body  undergoing 
tranaletlon  end  rotation,  it  haa  been  ahovn  (Ref.  D^that  thij  contribution  can  be  axpreeeed  in  terma  of  the 
inatantaneoua  rectilinear  and  angular  veloeitiea,  -v  and  u,  of  the  aolid.  Equation  (8)  then  redueea  to 
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■har.  K.  and  ara  raapaccivaljr  cha  rafiont  occupiad  by  cha  fluid  and  tha  aolid  body.  Tba  eoordinata 

tyacaa  tfanalataa  with  tha  aolid  but  doaa  net  rotate  with  it. 

vortical  regiona  at  any  finite  tiaa  level  after  tha  onset  of  tha  ■otlon  ara  of  finito  axtont  and 
a finit^  dlatancsa  frow  tha  solid  boundarioa.  Tha  region  of  integration  ia  tharafora  finite, 
a,  aa  r goes  to  infinity,  tha  integrands  in  Bq.  (9)  vanish.  Tha  fraastraaa  boundary  condition  ia 
'titfiad  at  infinity. 

Equation  (7),  togatbar  with  appropriats  boundary  conditiona,  can  ba  racaat  into  an  intagral  rapratan- 
, ion  for  tha  vaccor  potantial  T.  Tor  two-diaonaional  flowa,  T haa  only  ona  oon-vaniahing  coaponant  in 
tha  diraction  parpandicular  to  tho  piano  of  tha  flew,  Thia  cofoaant  la  faaillarly  intarpratad  aa  tha  atraaa 
function  T.  Aa  integral  rapraaantation  for  T ia  (Baf.  10): 
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Tha  aacoad  iaga,tal  ia  Bq.  (10)  agaia  rapraaanct  tha  eontributioa  of  tha  boundary  condition.  Tor  a 
aalld  hady  traaatatiat  ia  aa  iaflaita  flaid,  thla  aoatrihatioa  ia  aiaply  («sa  f).k,  whara  k ia  tha  unit 
wactor  ia  tha  diaactiaa  parpaadiealar  to  tha  floa.  Iqaatioa  (10)  kacoaaa  for  thla  caaa 


T (t.  t)  • r ( «(?  , t)  la  -i^l  a ♦ (^  a ?)•  t (11) 
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Bquatioaa  (10)  aad  (11)  ara  valid  ia  tha  phyaieal  plaaa  (a,  y)  and  in  any  plana,  aay  ((,n),  ralatad  to 
fha  phyaieal  plaaa  by  a eoaforaal  traaafanatlaa,  with  tha  atraaa  functiaa  T aa  invariant  uadar  tha 
traaaforaatioa.  Tha  vortieity  la  tba  (a,y)plana  ia  ralatad  to  that  ia  tba  ((,n)  plaaa  by  tha  (ollowinB 
forwula 

M(C,n  . t)  • hVm,  y,  t)  (12) 

:diora  ■ la  tha  aeala  factor  for  tho  traaaforaatioa  aad  ia  dafiaad  by 
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Bqwatioa  (11)  ar,  aqwivalaatly,  Bq.  (t)  foiaa  tha  haaia  of  aarliar  atudiaa  of  tha  iataira- 

diffaraatial  apTraath  hy  tha  yraaaat  aathara  aad  hy  athart.  It  hat  baaa  poiatad  oat  (Baf.  •)  that  tha 
aparatiaa  aoaat  af  aaaariaal  praaadaraa  haaad  aa  thaaa  aqaatiaaa  eaa  ha  high,  far  asaapla,  aappaaa  tha 
vortiaal  rapiaa  af  dM  flaw  caataiaa  ■ irid  pwiata.  Iqwatioa  (11)  la  rapraaaatad  hy  tha  alpahraie 
aqaatiaaa 

\(t)  a jli  ♦ •4.  i-  I,  (U) 

^ra  i rafara  ta  a grid  paiat  at  ahiah  tha  walwa  af  f ia  ta  ha  taapatad,  J rafara  ta  a hauadary  paiat  or 
a vortiaal  paiat,  aao  •MMtrit  aoaffitUatt  qfwaa  walooa  dapaad  aa  tha  ralativa  peaitioa  af  tha  peiata 
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i «nd  j|  «nd  G.  are  constant  values  corresponding  to  the  last  ter*  of  Eq.  (11).  The  computation  of  each 
stream  function  value  using  Eq.  (14)  requires  about  2N  algebraic  operations  (N  multiplicatons  and  N 
additions),  provided  that  the  geometric  coefficients  f..  are  known  and  need  not  be  computed..  The  use  of 
Eq.  (14)  to  evaluate  stream  function  values  at  all  the  grid  points  requires  a total  of  2N^  algebraic 
operations.  It  is  well  known  that  an  efficient  finite*'dif ference  or  finite-'clement  solution  of  the  two- 
dimensional  form  o^  Eq.  (7)  for  a field  containing  N grid  points  requires  a much  soMller  number  of 
operations  than  2N  . The  use  of  the  integral  representation  in  the  form  of  Eq.  (11),  or,  equivalently, 
Eq.  (9),  is  therefore  computationally  efficient  only  if  the  confinement  of  the  solution  field  to  the  vortical 
regions  results  in  a substantial  reduction  in  the  number  of  grid  points  required. 

For  many  high  Reynolds  number  flows,  the  vortical  region  constitutes  only  a small  portion  of  Che 
entire  flowfield  so  that  a substantial  reduction  in  Che  needed  number  of  grid  points  is  indeed  accomplished 
by  the  use  of  Eq.  (11).  It  should  be  emphaaited,  however,  that  even  if  the  reduction  in  the  number  of  grid 
points  is  not  substantial,  the  incegro-dif farential  approach  still  offers  highly  efficient  numerical 
procedures . 


One  method  of  accomplishing  superior  computational  efficiency  is  to  use  Eq.  (11)  only  in  the 
computation  of  stream  function  (or  velocity)  values  at  grid  points  surrounding  the  vortical  regions.  Once 
this  is  completed,  the  computation  of  stream  function  values  at  grid  points  in  the  vortical  regions  can  be 
carried  out  using  a finite-difference  or  a finite-element  method.  In  this  sMnner,  the  superiority  of  the 
integro-dif ferential  approach  is  assured  since  an  optianns  method,  including  a method  that  may  appear  in  the 
future,  can  be  incorporated  into  the  numerical  procedure  and  the  solution  field  can  be  confined  to  the 
vortical  regions. 

Another  method  of  improving  the  computational  efficiency  is  to  utilise  a flowfield  segmentation 
technique  (Ref.  11).  This  technique  has  been  suggested  in  conjunction  with  the  velocity  integrals, 
Eqs.  (8),  and  (9)  and  tested  on  the  problem  of  a flow  past  a finite  flat  plate.  The  application  of  this 
technique  in  conjunction  with  the  stream  function  formulation  is  outlined  here. 

The  flowfield  is  segmented  into  compartments  each  containing  a suitable  number  of  grid  points.  With 
known  vorticity  values  at  all  grid  points  at  a new  time  level,  new  stream  function  values  at  grid  points 
located  on  the  boundary  of  the  compartments  are  computed  from  Eq.  (14).  Values  of  the  nonaal  gradient  of  the 
stream  function  at  the  compartment  boundary  are  computed  from  an  equation  similar  to  Eq.  (14)  derived  by 
differentiating  Eq.  (11)  and  representing  the  result  by  algebraic  equations. 


Stream  function  values  at  grid  points  located  interior  of  the  compartments  are  now  computed  using 
Eq.  (10).  Suppose  the  vortical  region  is  segmented  into  several  compartments  each  containing  P number  of 
grid  points,  with  Q number  of  grid  points  located  on  the  compartment  boundary.  Equation  (10)  is  now 
represented  for  each  compartment  by  the  algebraic  equations: 
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where  the  first  suMStion  is  over  sll  grid  points  of  the  compartment  and  the  second  suosation  is  over  all 
grid  points  on  the  compartswnt  boundary.  Tha  use  of  Eq.  (IS)  to  compute  each  atream  function  value  requires 
2(P^2Q)  algebraic  operations.  This  operation  count  ia  clearly  substantially  smaller  than  2N  if  N is 
large.  If  the  solution  field  is  segmented  into  a larger  number  of  compertments,  then  the  operation  count 
2(P^2Q)  dacreaaas.  On  the  other  hand,  with  an  increasing  number  of  compartments,  the  total  number  of  grid 
points  located  interior  of  all  the  compartmenta  dacreaaea.  Aa  a result,  there  are  fewer  grid  points  in  the 
solution  field  for  tdiich  the  segmentation  technique  offers  a computational  advantage.  At  the  same  time,  the 
number  of  grid  points  located  on  compartment  boundaries  increaaea  with  the  number  of  compartments;  and  on 
thaaa  grid  points  not  only  the  stream  function  value  but  also  its  normal  gradient  must  be  computed.  There 
exists,  therefore,  an  optimum  number  of  compartments  for  each  problem  considered.  This  optimum  number 
depends  on  tha  total  number  of  grid  points  in  the  vortical  region  and  the  geometric  distribution  of  these 
grid  points.  It  ia  not  difficult  to  show,  however,  that  the  optimum  number  of  compartments  in  general 
increases  with  the  total  number  of  grid  points.  Consequently,  the  factor  of  reduction  in  cha  operation  count 
also  incraaaaa  with  tha  total  number  of  grid  points.  For  axampla,  with  1089  grid  points  arranged  in  a 
33  X 33  array,  the  operation  count  can  be  reduced  by  a factor  of  about  2.3  with  the  segmentation  technique. 
With  4225  grid  points  arranged  in  a 65  x 65  array,  the  operation  count  can  be  reduced  by  a factor  of  about 
4.6. 


Substantial  reductions  in  cha  total  operation  count  for  problems  requiring  a large  number  of  grid 
points  ia  a highly  significant  feature  of  the  flowfield  segmentation  technique  since  it  ia  praciaaly  for 
thaaa  problama  that  axcaaaiva  counter  time  prasanta  a aarioua  obatacla.  An  additional  highly  significant 
feature  of  this  technique  is  tha  inharant  flexibility  it  offers.  Tha  compartmenta  utilised  in  cha  solution 
procedure  can  be  of  any  arbitrary  site  or  shape.  Once  the  atream  function  (or  velocity)  values  are 
aatabliahad  on  tha  boundary  of  a compartment,  tha  cooputacion  within  that  compartment  can  be  performed 
independently  of  that  in  ocher  comparCmanta.  Conaaqutntly , cha  flowfield  aagmanCaCion  technique  is  well 
suited  for  parallel  programming. 

Sines  tha  comparCmanta  can  ba  choaan  to  ba  of  standard  aimpta  ehapae,  such  aa  ractaoglaa,  except  near 
solid  boufidariaa  of  irragular  shapa,  various  affleiant  mathoda  can  ba  adopted  for  tha  interior  point 
eo^utation.  That  it,  on#  does  not  always  nasd  to  uss  Iq.  (IS)  for  tha  intarior  point  computation.  Rathar, 
highly  affleiant  mathoda  such  aa  tha  faat  Fourier  transform  mathod  that  are  particularly  well  suited  for 
aimpla  boundary  gaoawtriaa  can  be  utilised.  If  Eq.  (15)  ia  used,  chan  the  number  of  geometric  coefficients 
that  antera  tha  computation  proeadura  ia  drastically  reduced  by  the  use  of  compartmenta  of  stands^  shapes. 
Rote  that  with  Eq.  (14),  the  number  of  geometric  coefficients  entering  the  computation  ia  N^.  Thaaa 
gaomatric  coaffielanta  are  Lndapendant  of  time.  Howavar,  If  M ia  vary  large,  then  tha  storage  of  tha  fr 
coefficianta  may  praaant  a problem,  and  tha  coaffielanta  may  need  to  ba  computed  rapaatadly  for  diffarant 
time  lavala.  With  gq.  (15)  and  standard  shape  for  compartments,  chs  numbsr  of  geometric  coefficients  is 


P(P*2Q).  This  numbtr  c«n  b«  further  reduced  to  P^2Q  If  the  compertisente  are  rectangles  and  unifonnly  spaced 
grid  la  uaed  in  each  compartment.  These  numbera  are  generally  small  enough  so  that  all  geometric 
coefficients  can  be  stored. 


As  discussed  earlier,  the  length  scales  for  the  potential  and  vortical  regions  of  the  flow  are 
generally  vastly  different  in  magnitude.  The  use  of  the  integral  representation  for  the  kinematic  aspect  of 
the  problem  makes  it  possible  to  confine  the  solid  field  to  the  vortical  regions  only.  However,  in  the 
general  viscous  flow,  the  length  scale  within  the  vortical  ragions  varies  greatly.  In  the  boundary-layers , 
the  length  scale  is  represented  by  the  boundary  layer  thickness  in  the  direction  perpendicular  to  the  solid 
surface.  In  the  wake  and  recirculating  regions,  the  length  scale  is  represented  by  the  body  dimensions.  The 
segmentation  technique  offers  a convenient  means  of  separating  these  parts  of  vortical  regions  in  the 
solution  procedure. 


The  above  features  of  the  segmenta^io;.  technique  are  expected  to  be  highly  important  for  three- 
dimensional  flows  requiring  the  use  of  very  large  number  of  grid  points.  For  the  present  studies  of  c%ro- 
dimensional  flows,  only  a few  of  these  features  of  flovfield  segmentation  are  incorporated  into  the  solution 
procedure.  It  is  noteworthy  that  the  segmentation  technique  is  useful  also  in  the  solution  of  the  kinetic 
aspect  of  the  problem.  In  the  present  study  of  the  impulsively  started  airfoil  problem,  the  solution  field 
is  divided  into  compartsMnts  in  the  kinetic  part  as  well  as  the  kinematic  part  of  the  computation. 

3.4.  Extraneous  Vorticity  Boundary  Condition 

The  solution  of  Sq.  (4)  requires  a knowledge  of  the  vorticity  values  on  the  boundaries  of  the  fluid  as 
a function  of  time.  The  establishment  of  new  vorticity  values  on  solid  boundaries  is  a necessary  part  of  the 
overall  solution  procedure.  In  most  of  the  previous  investigations,  one-sided  difference  formulas  have  been 
used  to  estimate  the  vorticity  values  on  solid  boundaries  from  the  no-'slip  condition  and  the  computed  stream 
function  (or  velocity)  values  at  grid  points  near  the  boundary.  The  computer  time  required  for  such 
calculations  is  very  small  compared  to  that  required  for  the  calculation  of  field  values  of  vorticity.  It 
has  bean  found,  however,  that  one**aided  difference  formulas  of  first-order  accuracy  tend  to  yield  stable 
solutions,  while  second-order  accurate  formulas  tend  to  give  unstable  results  at  high  Reynolds  numbers.  In 
fact,  for  some  problems,  second-order  formulas,  even  when  stable,  gave  leas  accurate  solutions  than  did  the 
first-order  formulas  (Referance  19).  The  use  of  firat-order  formulas,  however,  restricts  the  overall 
accuracy  of  the  solution.  In  particular,  it  has  been  shown  (Reference  8)  that  first-order  formulas  do  not 
permit  the  effect  of  tangential  pressure  gradient  on  the  local  boundary  vorticity  generation  to  be  accounted 
for.  Thera  axisted,  therefore,  considerable  uncertainties  regarding  the  correct  one-sided  formula  to  use, 
and  ragarding  tha  limitations  of  each  formula.  These  difficulties  are  independent  of  the  computer  time 
needed  for  tha  component  problem  under  consideration^  and  are  discussed  mure  fully  in  Reference  8.  They  are 
referred  to  as  the  extraneous  boundary-condition  problem  since  the  proper  boundary  condition  for  the 
physical  problem  is  that  for  the  velocity  (or  the  atreem  function).  The  specification  of  the  vorticity  on 
Che  boundary  of  the  fluid  domain  overspecifies  the  problem.  These  difficulties  of  extraneous  boundary 
conditions,  however,  are  not  peculiar  to  the  formulation  of  the  problem  using  the  vorticity  as  a field 
variable.  If  Bq.  (1)  is  to  be  solved  nuamrically  instead  of  Bq.  (4),  Chen  it  is  necessary  to  know  the 
pressure  on  Che  solid  boundaries  in  addition  to  Che  velocity. 

The  extraneous  boundary-condition  problem  is  absent  in  the  boundary  layer  problem  since  the  pressure 
in  a direction  normal  to  the  boundary  (vortical)  layer  is  practically  constant.  Thus  the  pressure  on  the 
solid  boundaries  la  equal  to  that  at  the  outer  edge  of  Che  boundary  layer  where  it  is  determined  by  the  outer 
flow,  For  the  general  viscous  flow,  the  vortical  region  is  not  necessarily  chin  and  the  above  simplification 
may  not  be  valid 

In  Ref.  8,  it  is  shown  that  the  prescribed  velocity  boundary  condition  imposes  a kinematic 
restriction  on  the  boundary  vorticity  values.  The  integral  representation  for  the  velocity  vector  then 
pnrmits  an  accurate  determination  of  the  boundary  vorticity  value  through  the  solution  of  an  integral 
equation. 


Consider,  for  example,  a solid  undergoing  translation  only.  The  second  integral  in  Bq.  (9)  vanishes. 
If  Che  vorticity  diaCribuCion  is  known  interior  of  the  fluid  domain  R.*,  then  the  first  integral  in  Bq.  (9) 
can  be  written  as  a sum  of  an  integral  over  R.'  and  an  Integral  over  the  solid  boundary  B on  which  a layer 
of  vorticity  C exiata.  For  two-dimeneionei  flow,  one  then  hee,  aince  v • 0 in  the  aSlid, 
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With  J know  in  t it  datarainata.  Equation  (16)  thnrnfort  in  t Fradhola  intngral  tquition 
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la  praviout  ttudiat,  tha  inttBral  aquation  (17)  it  approaiaattd  by  a aat  of  tiaulttnaoua  alEobrtic 
aqvatioat  coataiaiaE  C valuta  at  boundary  pid  peintt  at  unknown.  Thn  tat  of  nquttlona  it  nolvad  uain. 


iterative  methods.  In  the  present  study,  closed  form  solutions  of  the  integral  equation  (17),  subject  to  the 
restriction  (16),  are  obtained.  These  closed  form  solutions  make  it  possible  to  evaluate  values  on  the 
solid  boundaries  explicitly,  point  by  point.  It  is  found  that  these  explicit  calculations  yield  highly 
accurate  values  of  vorticity  on  solid  boundaries.  The  derivations  of  the  closed  form  solutions  are  presented 
in  Ref.  I and'will  also  appear  in  future  articles. 

4.  IMPULSIVELY  STARTED  AIRFOIL 

The  integro-differential  formulation  in  terms  of  the  stream  function  and  vorticity  is  particularized, 
transformed,  and  applied  to  the  study  of  the  time-dependent  incompressible  viscous  flow  past  an  impulsively 
started  9X  thick  symmetric  Joukowski  airfoil  at  an  angle  of  attack  of  15  and  a Reynolds  number  of  1000.  This 
problem  has  been  studied  previously  by  Mehta  (Ref.  22)  %dJo  mapped  the  fluid  domain  into  the  interior  of  a 
unit  circle  and  used  a finite-difference  siethod  to  so'.ve  the  transformed  vorticity  transport  and  Poisson's 
equations.  In  the  present  study,  the  geometry  of  the  airfoil  is  identical  to  the  one  studied  by  Mehta.  The 
vorticity  transport  equation  and  the  integral  representation  are  re-e ..pressed  in  Mehta’s  transformed  plane 
and  Mehta's  grid  spacings  are  adopted  in  order  to  facilitate  a compar.jon  between  the  integro-differential 
results  of  the  present  study  and  the  finite-difference  results  of  Mehta. 

The  major  features  of  the  numerical  procedures  used  in  this  study  is  described  below  for  the  component 
problems  shown  in  Fig.  1. 

The  grid  system  is  obtained  through  a transformation.  This  system  is  boundary-fitted  and  expanding. 
The  shape  of  the  9t  symaietric  Joukowski  airfoil  is  given  in  Figure  2 in  the  complex  plane  z ■ x ♦ iy.  By 
using  the  conformal  transformation 
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where  < ■ C and  y are  real  constants,  the  airfoil  surface  transforms  into  a circle  and  the  exterior  of 

the  airfoil,  i.e.,  the  fluid  domain,  is  isapped  in  the  x plane  as  Che  interior  of  the  ircle. 

For  the  present  problem,  the  values  of  C and  y are  taken  to  be  0.92A1635  and  -0.05214  respectively. 
The  coordinates  and  the  flow  variables  are  to  be  non-dimensionalized  with  respect  to  the  radius  of  the  circle 
and  the  freestream  velocity.  The  non-dimensional  chord  length,  L,  of  the  airfoil  is,  according  to 
equation  (19),  3.71261277.  Ihe  airfoil  is  6.9996Z  chick. 

The  kinematic  part  of  the  computation  is  performed  in  the  x plane.  Anon-conformal  transformation  is 
utilized  to  obtain  an  efficient  distribution  of  grid  points  for  the  kinetic  part  of  Che  computation.  The 
coordinate  in  the  < plane  is  not  laodified.  The  r coordinate  is  "stretched"  by  the  relation  (Ref.  22) 

p ■ t«nh'^(1.996590918r  - 1.032563339)  + 2.8  (20) 

This  transformation  maps  the  annular  region  between  r ■ 0.02  and  r ■ 1 in  the  x“plane  into  the  interior  of 
a unit  circle  in  the  p-8  plane. 

The  grid  system  in  the  p-0  plane  is  formed  by  using  Ap  ■ 1/47  and  A0  “ n/40.  The  coordinates 

of  the  grid  point  i,j  are 

p ■ 1.0  - (j  ' l)Ap  and  0 -(i  - 1 ) A0 

Mote  that  j • 1 corresponds  to  the  airfoil  surface  and  j • 48  (p  “o)  corresponds  roughly  to  a circle  of 
radius  13  L in  the  physical  plane.  The  corresponding  grid  system  is  shown  in  Fig.  2. 

The  vorticity  distribution  ionediately  after  the  start  of  the  airfoil  motion  consists  only  of  a vortex 
sheet  on  Che  airfoil.  The  determination  of  the  initial  vorticity  distribution  does  not  require  a special 
procedure.  Rather,  the  procedure  outlined  inSection  3.4.  is  applicable.  With  the  vorticity  everywhere  zero 
in  R ' , the  vector  f,  defined  by  Eq.  (17),  is  simply  -2xv  . With  this  value  of  Eq.  (16)  can  be  solved 
numerically  to  give  the  initial  distribution  of  vorticity  (vortex  strength)  on  the  airfoil  surface. 
Alternatively,  Che  closed  form  solution  of  Eq.  (16),  subject  to  the  condition  (18),  is  recognized  to  be 

C - 2v^Sin(0*a)  (21) 

in  Che  transformed  plane  The  corresponding  vortex  strength  in  the  physical  plane  t is  equal  to  C /H,  H 
being  the  scale  factor  defined  by 


The  vorticity  transport  equation  expressed  in  the  working  plane  in  terms  of  non-dimensional  variables 
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wh«Tt  u ii  the  eorticity  in  thn  phyeienl  ptnne. 
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Th*  vortlcity  transport  aquation  (23)  it  approaiaata^  by  an  iaplicit  finita-difftranca  aquation  uainq 
a 3-point  backward  difftrancinp  to  tt^rtaant  tba  tiat  dariuativa: 
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Tht  coaffieianta  I,  D,  t and  I contain  darivativat  af  T and  art  tiaa-dapandant.  In  tht  pratant  work,  tha 
coafficianta  art  awaluattd  at  tha  tiat  laval  a. 


Iha  aalaetion  of  eantral  diffaraaeiaf  for  tba  apaea  darivativaa  of  vorticity  waa  aada  afttr  coaparinp 
coaputad  raaulta  for  tba  first  faw  tiat  atapa  otini  eaatral  diffaraneina  with  thoaa  using  upwind  diffaroncing 
and  using  Arakawa  diffaraneing  (Baf.  23).  Be  tigaificaat  daviationa  batwaan  tha  nuatrieal  valuta  of 
vorticity  was  rawaalad.  Mahta  (Baf.  22)  usad  Arakawa  diffaraneing,  which  raquiraa  aora  aritliatcical 
optrationa  than  eantral  diffaraneing,  to  tnsurt  total  vorticity  constrvation.  Sines  tha  praaant  aathod  of 
dattmining  aurfact  vorticity  distribution  conaarvtt  tha  total  vorticity,  tha  uaa  of  eantral  diffaraneing 
requiring  ftwar  arithaatie  oparaciona  is  tppropriaCa.  Saall  aaniinta  of  apurioua  valuta  of  vorticity  was 
found  to  occur  in  portions  of  the  outer  fiald  uhara  tha  vorticity  values  art  aspaetad  to  ba  nagligibla.  Such 
occuraneat  art  perhaps  attributabla  (Bat.  24,  23)  to  tha  larga  grid  spacing,  and  hanct  call  Btynolda 
nuabara,  uaad  in  tha  pratant  study  io  ragions  af  tha  flew  far  froa  the  airfoil  aurfaca.  Tba  ate  af  upwind 
diffaraneing  in  plaea  of  central  diffaraneing  raaavaa  this  aaeaaloua  solution  behavior.  Arakawa  diffar- 
ancing  waa  not  found  to  bo  affactiva  in  raaoving  this  tporioua  vorticity.  Sines  apurioua  vorticity  values 
are  vary  aoall  and  tinea  thay  do  not  aariousty  affect  tha  solution,  it  waa  decided  to  use  eontral 
diffaraneing  which  is  higher  orderad  than  upwind  diffaraneing. 

A nodifiad  siiuogly  inplicit  procadura  (SIF)  (Baf.  2t)  was  uaad  to  tolva  tha  tat  of  aiwultanaous 
algebraic  aquations  (24)  for  vorticity  valuat.  at  tha  tias  laval  nvl.  Tbit  aathed  is  talactad  awoogtt  several 
conpating  wathodt  on  tha  batia  of  its  aolutioa  affieianey  and  aeenraey.  Tba  ralativa  raquiraaantt  of 
cowputar  tiaa  using  the  various  nathoda  ware,  few  the  first  faw  tiat  steps,  (a)  Modified  SIF  - lOOt, 
(b)  SIP  - 220S,  (e)  Altamating  Diraetion  Iaplicit  lABX)  aatbod  - 1S4S,  (d)  Suecasaivo  Lina  Balanation  - 
2S4X,  and  (a)  Suceaativa  Foint  Balasation  Method  (naod  in  Bafaraaea22)  - SSSS.  Tba  AMX  atthad  it  faster 
than  tha  SXF  aathod;  but  did  net  yield  aceoptablo  aalation  naing  tba  saaa  tiat  atap  aiaaa  aaplayod  by  tba  SIF 
aathod. 


Tha  SIF  aathod  wot  proposed  by  Stoat  (Baf.  IS).  Lin  at  al.  aaplayod  tba  aathod  aaeeoasfully  to  obtain 
tolutiona  for  viteoua  float  past  a eyliadar.  Tha  baaie  oaacapta  af  this  oatbad  m applied  to  tba  pratant 
problaa  is  aa  followa.  Tba  eoafficiant  antrin  (or  Ig.  (14)  ia  aparao  aad  eantala  anly  fiva  naa-soro 
diagonals.  With  a direct  aliainatiao  praaadura,  oaa  waoXd  faetoriaa  tha  eaaffieiaat  antrin  into  tba  product 
of  a looor  aad  an  oppar  criaagular  aatria,  or  oao  a preeadors  ttat  ia  aaaaatially  oqaivalant.  In  such  eaaaa 
tha  triangular  ontricat  gaaeratad  are  not  ia  gaaaraX  aparaa.  With  the  SXF,  tha  eaaffieiaat  antrin  ia  altarad 
la  such  a aaaaar  that  the  altarad  antrin  is  aaaily  factanad  iata  a praduet  of  a Xoaar  and  aa  upper  triaa^ular 
aatria  and  that  the  nntriaoa  each  hovo  only  tbrao  anwaare  olanaata  ia  aaeh  raw. 


With  tha  SXF  aathod,  tha  aolutioa  field  ia  booadad  by  tha  solid  aarfoaa  (p  • 1)  and  a eoaataat  p line, 
■■y  F * Pi  idiieh  is  auffieiaatly  roaovod  froa  tha  ato-asgligibla  vorticity  rogioa.  Tba  valua  of  p , it 
gaoarally  datarainad  by  tba  ontaat  of  tba  woks  rogioa.  At  large  tiaa  lavsls,  tbs  wsks  antaads  far  bohiaa  tht 
airfoil  while  tbs  vortical  ragiont  to  tha  tidaa  and  ahead  of  tha  airfoil  aataad  only  a rolativaly  short 
dittaoea  froa  tha  airfoil.  Aa  a raault,  a larga  auahar  of  data  poiata  uhara  vorticity  valosa  ora  nagligibla 
antar  into  tha  solutian  procadura.  Tba  oedificatiaa  of  tha  SXF  aathod  latraduaad  ia  tho  praaant  study 
centiata  of  a division  of  tha  aolutian  fiald  into  caworal  ea^ortaanta  alaag  anna  tent  S lines,  tho  individual 
applicatiaa  of  SXF  aathod  in  each  caopartanat,  aad  tha  ueo  at  auccasaivo  liaa  ralaastioa  atthad  for  obtaining 
vorticity  valuac  aloH  the  dividing  coostoat  S liaas.  With  the  solution  field  dividad  ia  thia  onaasr,  aaeh 
solution  block  ia  banadaS  by  two  caaatsntS  liaaa,  a part  at  tha  p*  1 liaa,  aad  p ■ ft  liaat,  p.  can  ba 
different  for  difforant  blocha.  Thorsfaro  nuch  fswar  data  paiata  at  whiah  vorticity  watuoa  ato  oHiisi^i* 
antar  iata  tha  ecoputatioa.  In  tha  praaant  atudy,  tha  aalntion  fiald  ia  dividad  into  four  caopartaanta.  Tba 
location  of  tha  caostaat  • dividiat  liaaa  aad  p.  vaXoaa  sro  bapt  flanibla  sad  ara  ahoagad  in  asaardaaca  with 
tha  ahopo  of  tha  aaa-nagligibla  vorticity  ragiah  as  tiaa  pcogrosaas. 


Tba  vortsn  atrangths  an  tba  airfoil  aurfaca  ara  caapotod  using  oa  anplicit  foranla  in  tha  traaafemad 
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Tbs  value  of  tba  airfoil-bouadary  vavtlaity  ia  tho  yhpaiaal  pleas  at  a grid  point  located  oa  the 
airfoil  bauadary  is  tlMa  abtaiaad  by  dividing  ( by  Who  aealo  factor  ■ aad  by  the  physiaal  distoaca 
carraapsadlag  ta  >iSp  nsorost  tha  airfoil  aorfaoa.  The  oorfaaa  vorticity  values  ara  ro-aslsnlatad  for  each 
itarstioa.  Maw  baandary  vortiaity  vatoas  Car  Mm  aobaagaant  itarotlaa  arc  abuiaad  aaiag  a ralaaatiaa 
praaadura  and  nsod  in  tha  IXF  praaadara  for  tha  aalialatSaa  af  vaeticity  valaaa  at  grid  poiata  ooay  froa  tha 
baandary. 
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The  coaputation  of  the  atreaa  function  ia  accoapliahed  uaing  the  flovfield  oegaentation  technique 
deacribed  in  Section  3.3.  The  flovfield  ia  divided  into  aevcral  coapartaenta.  The  boundary  of  the 
coapartaenta  are  the  0*0  line  and  two  p ■ conatant  linea.  Streaa  function  valuea  on  conatant  p 
coapartaent  boundariea  are  coaputed  uaing  Eq.  (14)  in  the  tranaforaed  plane.  Streaa  function  valuea  at  the 
interior  pointa  of  each  coapartaent  are  then  calculated  uaing  the  atrongly  iapHcit  owchod  deacribed 
earlier. 


Inatantaneoua  preaaure  diatributiona  and  ahear  atreaa  on  the  airfoil  aurface  are  deterained  froa  Che 
coaputed  vorticity  valuea  and  gradienta  on  Che  aurface  in  accordance  wich  che  following  equation# 

If  - p'*  1^ 

and 
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Note  that  the  direction  of  the  normal  coordinete  ia  away  froa  che  fluid  doa^ln.  Equation  (26)  ia  a 
apecialisation  of  Eq.  (1)  to  Che  aolid  boundary  S where  the  "no-alip'*  condition  appliea.  The  preaaure  and 
akin  friction  coefficienta  C and  Cf  are,  froa  Bqa.  (36)  and  (37),  eapreaaible  aa  a function  of  6 aa 
follow*:  P 

p*l  o ^ p*l 

and  ^ 

C.(e)  - iru*|  (29) 

***  p-1 

The  force  and  aoaent  are  readily  obtained  froa  the  coaputed  valuea  of  C ^)  and  C^(q),  by  nuaerical 
quadrature.  In  the  preaent  work,  the  valuea  of  / Sp  are  obtained  by  nuaeri<fal  differentiation  uaing  che 
valuea  of  ui^  ^ and  2* 

Nuaerical  reaulta  cdtaioed  for  the  airfoil  problem  are  preaented  in  Figa.  4 to  9.  The  chordwiae 
diatancea  froa  the  leading  edge,  and  che  load  coefficienta  C,  , and  ahovn  in  che  figure#  are 
noraaliaed  wich  reapect  to  the  chord  length.  The  Icynolda  nuaber  of  lOOCTia  baaea  on  the  chord  length  and  the 
freeatreaa  velocity.  All  other  quantitiea  ahown  were  noraaliaed  wich  reapect  to  che  freeatreaa  velocity  and 
the  radiua  of  the  unit  circle.  The  reference  time  ia  the  radiua  of  the  circle  divided  by  Che  freeatreaa 
velocity. 

A check  on  the  program  coding  waa  carried  out  by  obtaining  the  aolution  for  a flow  paaC  a circular 
cylinder  at  a Kcynolda  nuaber  of  40.  The  aolutiona  were  carried  out  up  to  a diaenaionleaa  time  level  of 
t * 9.5,  the  reference  time  being  the  cylinder  radiua  divided  by  the  freeatreaa  velocity.  The  drag 
coefficient  at  Chia  time  level  ia  found  to  be  1.67  which  ia  in  reaaonable  agreement  wich  an  experimental 
value  of  1.5  for  ateady  atate  flow.  The  coaputed  preaaure  diatribution  on  che  cylinder  ia  coapared  with  the 
experimental  valuea  (Bef.  2S)  in  Figure  3. 

Additionally,  for  the  airfoil  problem,  the  acceptibility  of  block  aubdiviaion  waa  checked  out  by 
comparing  Che  nuaerical  reaulta  for  a few  firat  time  atepa  of  che  following  caaea:  (1)  Succeaaive  Point  and 
Line  Relaxation  Procedure#  without  block  aubdiviaion#  (2)  Strongly  Implicit  Procedure  (SIP)  with  no  flow 
field  aubdiviaion  (3)  81F  method  with  four  block  aubdiviaion#  and  for  three  different  locationa  of  block 
boundariea  for  the  block  encloaing  the  trailing  edge.  The  excellent  agreeaenc  of  computed  vorticity  valuea 
for  all  theae  caaea  indicated  the  acceptability  of  block  aubdiviaion  technique.  The  major  aapecta  of  the 
flow  f^enoaena  paaC  che  iapulaively  acarted  airfoil  at  an  angle  of  attack  of  15  and  at  a Reynold#  number  of 
1000  arc  deacribed  below. 

At  Ciaa  t * 0,  Che  flow  around  Che  airfoil  correaponda  to  that  of  a non*-circulaCory  potential  flow. 
After  the  i^ulaive  atart  the  rear  atagnation  point  (R8F)  aovea  very  rapidly  to  the  trailing  edge  (TE)  froa 
it#  potential  flow  location  of  X.  * 0.9808.  Thia  rapid  moveaent  of  che  rear  atagnation  point  ia  accompanied 
by  a aaall  aeparation  bubble  locAed  between  tSP  and  T?  at  Ciaa  level#  0.003  and  0.004.  After  the  RSP  haa 
reached  the  TB,  plot#  of  aqui*vorCicicy  contour#  ahov  **cttrling  up"  of  vorticity  field  near  trailing  edge  with 
poaitive  vorticity  valuea  in  Che  curled  up  region#.  The  curling  up  of  vorticity  field  ia  indicative  of 
vortex  roll-up  and  che  formation  of  atarting  vortex.  The  vorticity  in  thia  atarting  vortex  diffuaea  rapidly. 

Aa  adverae  preaaure  gradient  on  the  upper  aurface  becoaaa  noticeable  at  tiae  C * 0.068.  At  time 
t * 0.404  thia  adverae  preaaure  gradient  ia  aore  pronounced  (Fig.  4)  and  chia  lead#  to  the  birth  and  growth 
of  the  firat  aeparation  bubble  with  ctockwiae  flow.  Thia  bubble  eventually  extend#  over  alaoac  che  entire 
upper  aurface  of  che  airfoil,  at  tiaa  t ■ 6.228  (Fig.  5a).  At  the  very  next  tine  level,  t • 6.74,  ita  re- 
attachment  point  aeparacea  from  the  airfoil.  In  Figure  5,  atreaalinea  are  ahown  around  che  airfoil  for 
aeveral  tiaa  level#.  The  valuea  of  the  atreaa  function  are  in  the  range  -0.48  to  0.  6 in  atepa  of  0.04.  Aa 
the  bubble  aiae  iocrcaeee  the  preaeure  gradient  on  the  upper  aurface  decreaaea  aa  aeea  froa  Fig.  4.  The  lift 
coefficient  which  haa  been  rapidly  decreaaing  after  the  iapulaive  atart  (Fig.  6),  begina  to  increaae  wich  the 
growth  of  Che  eleckwiae  bubble  A,  aa  aeen  froa  Fig.  7.  Thia  increaae  in  lift  ia  aaaociated  with  che  fact  that 
the  attaehed  bubble  effectively  increaae#  the  camber  of  the  airfoil.  After  the  burating  of  tha  bubble  at 
time  t * 6.74,  the  lift  coefficient  increaaea  at  much  alower  rate  and  it  atarta  falling  off  at  tiaa 
t • 9.0  (Fig.  8).  The  burating  of  the  bubble  eauaea  the  drag  coefficient  to  increaae  (Fig.  8),  aince  che 
wake  width  ia  increaaed  for  che  aeparated  flow. 

Aa  adverae  preaaure  gradient  JSFthejrevereed^^nM  on  the  upper  aurface  from  X 0.82  to  X 0.6  ia 
noticeable  in  Fig.9.A  countercloekwiae  Uie^ie,  h,  appear#  ia  thia  region  and  growa  witif  time.  The  |rowth  of 
thia  countercloekwiae  bubble  ceuaea  the  lift  coefficient  to  fall  off  rapidly  (Fig.  8).  Froa  ciaa  t • 6.74 
to  t • 17.62  there  ie  no  rear  atagnation  atreaalioe  and  the  rear  atagnation  point  becoaaa  che  aeparation 
point  of  a countercloekwiae  bubble,  C (Fig.  56).  The  growth  of  bubble  8 partition#  the  recirculating  region 
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A into  regions  A'  and  D (Pig.  3b).  Also  bubbles  B and  C merge  together  and  then  lift  off  from  the  surface  at 
time  t • 17.748  (Fig.  5c).  AC  chis  time  level  region  0 becomes  an  attached  clockwise  bubbla,  and 
subsequently  at  time  t * 18.26  a counterclockwise  bubble,  E,  appears  within  the  attached  bubble  D.  Growth 
of  this  counterclockwise  bubble  E divides  bubble  D into  two  parts  (Pig.  5d).  As  time  progresses,  bubble  D* 
of  Pig.  3d  moves  towards  the  trailing  edge  and  bursts.  After  the  bursting  of  D',  the  bubble  D grows  rapidly 
in  site  and  at  time  t • 23.428  it  extends  over  alsKist  the  entire  upper  surface  of  the  airfoil. 

The  present  results  are  in  good  agreement  with  those  of  Hehta  (Ref.  22)  until  the  bursting  of  the 
first  separation  bubble.  After  this  the  sequence  of  occurence  of  bubbles  and  chair  bursting  or  lift  off  from 
surface  are  different  between  these  two  studies.  The  key  differences  are  the  following:  (1)  In  the  present 
study,  whenever  there  is  no  rear  stagnation  streamline,  a small  counterclockwise  trailing  edge  bubble  is 
always  present  and  the  rear  stagnation  point  becomes  the  separation  point  of  this  bubble.  This  ia  not  found 
true  in  Hehta*a  study.  (2)  After  bubbles  B and  C merge  they  stay  merged  and  lift  off  from  the  surface 
together.  In  Mehta'e  study  they  separate  again  and  only  bubble  C lifts  off  Co  surface.  Bubble  B stays  on  the 
surface  to  play  Che  roll  of  bubble  E of  the  present  case. 

The  following  general  concluaions  are  similar  between  the  two  studies.  (1)  iMsdiately  after 
impulsive  start  the  rear  stagnation  point  moves  very  rapidly  to  the  trailing  edge.  (2)  Bubsequantly  a 
"starting  vortex"  is  visualised  through  concentric  equivorticity  lines.  (3)  The  formation  of  a separation 
bubble  is  preceded  by  an  adverse  pressure  gradient  (4)  Clockwise  bubbles  extend  or  move  towards  trailing 
edge  and  burst  (5)  Anticlockwise  bubbles  either  lift  off  the  surface  or  open  up  to  streamlines  from  above 
Che  surface  (6)  The  lift  increases  with  increase  in  the  site  of  attached  clockwise  bubbles  and  dacraases 
when  attached  anticlockwise  bubbles  grow. 

3.  OSCILLATING  AIRFOIL 

The  inCegro-dif ferencial  formulation  in  terms  of  Che  velocity  vector  and  the  vorcicity  is  utilitad  in 
the  study  of  incompressible  flows  past  a 12X  Chick  symmetric  Joukowski  airfoil  oscillating  in  pitch  about  a 
pitching  axis  located  1/4  chord  from  Che  leading  edge.  The  major  features  of  Che  numerical  procedures  are 
described  below  in  terms  of  the  component  problems  shown  in  Fig.  1. 

The  solution  field,  already  confined  to  the  vortical  region,  is  divided  into  an  inner  region  and  an 
outer  region.  In  Che  inner  region,  the  grid  system  consists  of  grid  points,  or  nodes,  that  are  not  uniformly 
spaced.  The  nodes  represent  vortices  of  triangular  elements.  In  Pig.  10  is  shown  the  specific  grid  system 
for  Che  inner  region  used  with  the  12Z  thick  airfoil.  This  system  has  232  nodes,  including  48  nodes  on  the 
airfoil  surface  and  62  nodes  on  the  outer  boundary  of  Che  inner  region*  This  system  gives  354  triangular 
elements.  The  grid  system  for  the  outer  region  consists  of  uniformly  spaced  grid  points  forming  a 
rectangular  array.  The  outer  region  overlaps  the  inner  region  by  one  layer  of  elcmenta.  That  ia,  in  the 
outermost  layer  of  triangular  elements,  each  pair  of  elements  form  a rectangle  idiose  boundaries  coincide  with 
the  rectangular  grid  lines  of  the  outer  region.  The  grid  spacings  usad  for  the  airfoil  are  Ax  ” 0.05  and 
Ay  ■ 0.023.  This  grid  system  for  the  outer  region  is  not  shown  in  Fig.  10.  The  number  of  grid  points  in  the 
outer  region  is  large.  In  Che  oscillating  airfoil  case,  w»re  Chan  3500  points  are  used  in  the  outer  region. 

A f inite-element  method  is  used  to  treat  the  kinetics  of  the  problem  in  the  inner  region.  In  the 
outer  region,  fioite*dif ference  methods  arc  usad.  This  hybrid  finite  element*finlte  difference  method 
permits  the  grid  system  Co  "fit"  Che  airfoil  boundary  geometry  and  to  have  relatively  closely  spaced  grid 
points  near  Che  airfoil  boundary.  By  using  a relatively  smell  number  of  nodes  in  Che  inner  region,  the  siae 
of  the  relatively  complex  coefficient  matrices  involved  in  the  finite-elemant  method  ia  kept  small.  Moat  of 
the  solution  field  is  covered  by  the  f inice~dif ference  grid  which  leads  to  relatively  simple  coefficient 
matrices. 


Associated  with  each  node  "i"  in  Che  inner  region  at  (i,y  ),  there  is  a composite  linear  interpolation 
function  N^(x  ,y  ) with  the  property 


(30) 


where  (.  . ia  Che  Kronecker  delta.  The  interpolation  function  is  continuous  across  Che  element  boundaries 
(It  ia  tero  in  elements  not  containing  the  node  i)  chough  its  first  derivatives  are  not  necessarily 
continuous.  The  velocity  and  Che  vorcicity  fields  are  assumed  to  vary  linearly  within  each  element*  Thus 
one  writes 


u)  - (h.  N.  H.) 

).  k' 


(31) 


for  ch*  triangular  alaaant  with  noda*  i,  j and  k.  Uain,  Calarkin'a  procadura,  an  appreaiaation  to  tha  two- 
diaanaional  vortieitp  tranaport  aquation  ia 


"i^*’  y^^***)^  “ ® <32) 

Th*  abov*  waighting  procaaa  ia  appliad  at  all  noda*  in  tha  innar  ragion  aacapt  thoaa  on  tha  boundary  of  that 
ragion.  ly  tha  uaa  of  th*  divarganca  thaoran,  on*  obtain* 

v/  (V*iii)N.da  dy  • - V / ?H.  • ^ dxdy 
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The  line  integral  it  performed  over  the  boundary  of  the  inner  region.  Since  i it  not  a boundary  node,  thii 
line  integral  vaniahea  by  virtue  of  Eq.  (30).  By  the  uae  of  Eq.  (33),  the  integrand  in  Eq.  (32X  ia  re~ 
expreaaed  aa  one  containing  only  firat  derivativea  of  u) . 

Equation  (32)  can  be  written  aa  the  aum  of  component  integrala  over  all  alcmcnta  in  the  inner  region. 
Upon  introducing  Eq.  (31)  into  the  integrala,  each  of  the  component  integrala  in  the  aum  can  be  integrated, 
yielding  a ayatem  of  aimultaneoua  differential  equationa  of  the  form 

(H)  * (C)j.-]  ♦ (d)[u,|  - |e}  (34) 
Uaing  a backward  difference  acheme  for  the  time  derivative,  one  obtaina 


The  terma  on  the  left  hand  aide  of  Eq.  (32)  repreaent  reapectively  the  local  time  rate  of  change  of 
vorticity,  the  convection  of  vorticity,  and  the  diffuaion  of  vorticity.  The  coefficient  matricea  (h)  and  (o) 
are  depdendent  only  on  the  node  coordinatea  and  are  independent  of  time,  they  do  not  need  to  be  evaluated 
repeatedly  for  different  time  levela.  The  matrix  (c  ) repreaenta  convective  proceaaea  and  ia  dependent  on 
time  through  the  time-varying  velocity  diatributton.  To  avoid  iterative  computation  of  the  velocity  valuea 
for  each  time  atep,  the  velocity  field  ia  computed  only  for  the  old  time  atep  "n*'.  The  column  matrix  (E)  on 
the  right  hand  aide  of  Eq.  (35)  ia  alao  time  dependent.  It  repreaenta  the  contributiona  of  the  boundary 
nodea  of  the  inner  region.  Uith  known  vorticity  valuea  on  the  boundary  nodea,  (B)  can  be  explicitly  computed. 
In  the  preaent  work,  the  boundary  of  the  inner  region  conaiata  of  the  aolid  aurface  and  the  boundary  located 
within  the  outer  region.  (Recall  that  the  two  regiona  overlap).  Vorticity  valuea  at  boundary  nodea  that  are 
located  inaide  the  outer  region  are  computed  aa  a part  of  the  finite-difference  procedure.  Vorticity  valuea 
at  the  aolid  aurface  are  determined  aa  a part  of  the  kinematic  computation  aa  deacribed  in  Section  3.4. 

Since  the  matrix  (c  ) ia  tiom-^^pendent , the  coefficient  matrix  of  the  ayatem  of  aimultaneoua 
algebraic  equationa  for  the  unknown  a)  ” muat  be  inverted  at  each  time  atep.  It  ahould  be  noted,  however, 
that  fairly  large  time  intervala  can  be  uaed  in  the  aolution  procedure. 

In  the  outer  region,  the  vorticity  tranaport  equation  ia  treated  uaing  a two-point  backward  tiam 
differencing  acheme,  with  aucceaaive  point  overrelaxation  method  uaed  to  aolve  the  reaulting  algebraic 
equationa . 

The  velocity  values  are  computed  by  a numerical  quadrature  of  Eq.  (9).  The  firat  integral  in  Eq.  (9) 
ia  expreaaed  aa  a aum  of  component  integrala  over  individual  elementa.  (for  thia  purpoae,  each  rectangle 
bounded  by  grid  lines  in  the  outer  region  ia  treated  aa  a rectangular  element).  The  uae  of  interpolation 
functions  then  leads  to  exact  analytical  expreaaiona  for  the  component  integrala.  These  exact  expreaaions 
have  been  derived  for  interpolation  functions  of  all  degrees.  In  the  preaent  work,  linear  interpolation 
functions  are  uaed  for  the  triangular  eleawnta  of  the  inner  region  of  the  grid  ayatem.  The  vorticity  in  the 
rectangular  elements  of  the  outer  region,  however,  are  represented  by  e concentrated  vortex  filament  placed 
at  the  centroid  of  the  element  and  its  contribution  to  the  velocity  field  is  computed  in  accordance  with  the 
Biot-Savart  law.  The  evaluation  of  the  second  integral  in  Eq.  (9)  ia  simplified  by  considering  the  rotating 
coordinate  ayatem  and  the  nonrotating  coordinate  ayatem  to  coincide  with  one  another  at  the  time  t.  The 
second  integral  can  then  by  re-expreaaed  aa 


r • r 

2 J5(t)  X / 2 ISR  (J6) 
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Th«  integral  therefore  neede  to  be  evaluated  only  for  unit  angular  velocity  Q • 1 eince  it  ia  directly 
proportional  to  (i  . 

The  coaputetion  of  the  preaeure  and  akin  friction  coefficienta  are  carried  out  in  the  aenner  deacribed 
in  Section  (4). 

The  nuaerical  procedures  used  in  this  atudy  «ere  calibrated  by  treating  the  problea  of  a flow  past  an 
iapuleivaly  started  finite  flat  plate  at  a keynolde  nuaber  of  1000  and  at  laro  angle  of  attack.  Tha  results 
obtained  are  in  aacallant  agreaawnt  with  aarliar  raaulta  (laf.  $)  obtainad  uaing  finita-ditfaranca  aathoda 
inttaad  of  tha  praaant  hybrid  aathod  in  trtating  tha  vorticity  tranaport  aquation.  Tha  aolution  waa  carried 
to  a diwanaionlaat  tiaa  laval  of  2.4,  tha  rafaranca  tiaa  Isval  being  tha  plate  length  divided  by  tha 
fraaatraaa  velocity.  At  thia  tiaa  level,  tha  solution  varlaa  vary  slowly  with  tiaa  in  tha  vicinity  of  tha 
plats.  Tha  velocity  and  vorticity  profilaa  at  aidplata  obtainad  for  this  large  tiaa  lava!  era,  aa  aspacead, 
in  axcallant  agraaaant  with  tha  llaaiua  profile. 

At  an  additional  calibration,  tha  flow  paac  an  iapulaivaly  startad  circular  cylindar  at  a ksynolda 
nuabar  baaad  on  tha  cylindar  diaaater  of  40  waa  treated.  Tha  aolution  waa  carried  to  a diaanaionlaaa  tiaa 
laval  of  20.4.  At  thia  tiaa  laval,  tha  coaputad  drag  coefficient  variaa  only  in  tha  fourth  aignificanC  digit 
batwaan  aucesaaiva  tiaa  lavala,  the  tiaa  step  being  0.123.  Tha  coaputad  ttreaaline  and  conatant  vorticity 
contours  era  shown  in  Fig.  11,  in  tha  upper  and  lower  halves  of  tha  figure  rsapsetivaly,  for  this  large  tiaa 
laval.  hha  coaputad  praaaura  distribution  around  tha  cylindar  la  shown  in  Figure  3,  eoapirad  with  tha 
results  di-acribad  in  Section  4 and  axparlaantal  results  (taf.  21). 

Tha  caaa  where  tha  airfoil  ia  tat  into  translational  aotion  iapulaivaly  at  an  angle  of  attack  of  3** 
and  a Ksynolda  nwkbar  of  1000  waa  studied.  Tha  transient  solution  after  tha  onaac  of  tha  aotion  ia  carried  to 
a diaanaionlaaa  tiaa  laval  of  2.07,  tha  rafaranca  tiaa  being  tha  chord  length  divided  by  the  fraaatraaa 
velocity.  Tha  ttreaalina  pattern  at  a tiaa  level  of  0.S7  it  shown  in  Fig.  12.  Tha  adwarsa  praaaura  gradient 
on  the  upper  aurface  waa  not  sufficiently  strong  to  esuaa  flow  aaparation.  Tha  scrsaalina  pactarn,  however. 


i 
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shows  s thick  boundary  lsy«r  region  on  the  upper  surface. 

The  loads  on  Che  airfoil  are  very  high  iiMediaCely  after  the  impulsive  start.  They  decrease  rapidly 
at  small  time  levels,  up  to  a dimensionless  time  level  of  0.35.  Prior  to  this  time  level,  the  loads  are  due 
primarily  to  the  impulsive  start  effects.  Subsequent  to  this  tisM  level,  viscous  effects,  particularly  the 
thickening  of  the  boundary  layer,  results  in  further  changes  of  the  loads.  The  surface  pressure  decreases 
relative  to  the  trailing  edge  pressure.  The  upper  surface  pressure,  however,  decreases  at  a more  rapid  rate 
so  chat  the  lift  coefficient  increases  with  time.  Figure  13  shows  the  variation  of  the  lift  coefficient  C, 
and  of  the  drag  coefficient  C.  over  Che  time  interval  form  0.3  to  2.07.  Both  Che  lift  and  Che  drag 
coefficients  vary  relatively  slowly  at  the  time  level  2 07.  y 

A flow  past  the  airfoil  oscillating  at  a moderate  reduced  frequency  of  0.3  and  a teynolds  number  of 

1000  was  treated.  The  airfoil  is  set  into  translational  motion  impulsively  at  t ■ 0.  The  angle  of  attack  a 

is  given  by 

a - 3°  ♦ l”  Sin(0.6t)  (37) 

The  solution  was  carried  to  a non*‘dimensional  time  level  of  10.4.  The  airfoil  undergoes  a complete 
cycle  of  oscillation  during  this  period. 

The  main  feature  of  the  load  history  is  that  the  lift  and  moment  agree  qualitatively  with  the  time* 
dependent  linearised  potential  flow  solution.  The  calculated  lift  variation  is  in  phase  with  the  potential 
flow  result.  The  magnitude  of  the  lift  coefficient  is  lower  Chao  Che  lift  predicted  by  the  potential  flow 
analysis.  This  is  expected  since  at  the  moderate  leynolda  number  of  1000,  the  boundary  layer  on  the  upper 
surface  is  thick  and  its  displacement  thickness  results  in  a significant  change  in  the  effective  ahape  of  the 
airfoil . 

The  viscous  drag  is  much  larger  than  the  pressure  drag  for  the  present  case  of  low  angle  of  attack. 

This  viscous  drag  is  relatively  insensitive  to  the  pitching  motion  of  the  airfoil.  lumerical  results 

indicate  a very  small  separation  bubble  appearing  near  the  trailing  edge  on  the  upper  surface  during  the  down 
stroke  (not  at  maximum  angle  of  attack).  However,  the  site  of  this  bubble  is  very  small  and  the  presence  of 
the  bubble  does  not  affect  the  integrated  load  significantly. 

The  solution  obtained  at  the  end  of  the  moderate-frequency  oscillation  cycle  was  used  as  the  initial 
solution  for  a high  frequency  oscillation  computation.  The  flow  teynolds  number  is  1000  and  the  reduced 
frequency  is  3 for  this  case.  The  angle  of  attack  is  given  by 

a - 3°  ♦ 1°  Sin(6t)  (38) 

The  coapuCation  wa*  perforaad  for  aorc  than  two  cyclaa  of  oacillation.  Tha  tranaiant  affaet,  8ua  Co 
cht  initial  condition  uaad,  ia  praaant  during  chn  firac  cycla  but  bacaaa  vary  aaall  durlag  tha  aaeoad  cycla. 

The  load  hiatory  ia  praaancad  in  Figuraa  14,  IS,  16  and  17. 
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